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Abstract 

We construct new traveling wave solutions of moving kink type for a modified, 
driven, dynamic Frenkel-Kontorova model, representing dislocation motion under 
stress. Formal solutions known so far are inadmissible for velocities below a thresh- 
old value. The new solutions fill the gap left by this loss of admissibility. Analytical 
and numerical evidence is presented for their existence; however, dynamic simula- 
tions suggest that they are probably unstable. 

1 Introduction 

The first study of dislocation motion through a discrete model was performed by Atkinson 
and Cabrera pQ , followed by Celli and Flytzanis [4] and Ishioka [8 J . Atkinson and Cabrera 
PQ utilize a variant of the discrete sine-Gordon equation, which arises in the dynamic 
version of the Frenkel-Kontorova model. They seek traveling wave solutions corresponding 
to uniformly moving dislocations under applied stress, represented by a constant forcing 
term. To the present day, analytical progress on the forced problem has been limited, 
while some numerical conclusions have been drawn by Peyrard and Kruskal [TT]. In PQ, the 
trigonometric potential of the Frenkel-Kontorova model, responsible for the nonlinearity 
in the sine-Gordon equation, is replaced by a piecewise smooth potential with quadratic 
wells. A similar choice is made in jH [8]. The resulting equation for traveling waves 
reduces to a linear one, provided the solution satisfies an admissibility condition. This 
requires that to the right of a single transition point (the dislocation core), the solution 
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take values entirely in one of two quadratic valleys of the potential, and in the other 
valley to the left. Solutions of the reduced linear equation are found semi- analytically. 
Traveling waves are found to exist only if the stress and speed satisfy an algebraic relation, 
the kinetic relation of the dislocation. A remarkable prediction of this relation was that 
dislocations exceed the speed of shear waves at sufficiently high stress. This was recently 
confirmed experimentally [ID] . Another feature is the presence of multiple singularities 
and discontinuities in the graph of stress versus velocity, located at a sequence of special 
resonance velocities that accumulate at zero; see Fig. [TJ However, as emphasized by 
Earmme and Weiner [5] (see also [9]), below a threshold velocity, solutions of the linear 
problem violate the admissibility condition; this issue was not fully recognized in pQ. This 
rules out solutions in the entire interval containing singularities as inadmissible, since all 
resonance velocities are below the threshold velocity. 

The problem of existence of traveling waves at velocities below the threshold value has 
remained open to the best of our knowledge; this also applies to the results of [3j. This is 
the main issue addressed in the present paper. 

We present compelling analytical and numerical evidence that for velocities below the 
threshold value, there is a new type of traveling wave solution £ = x — Vt, which 
actually violates the usual admissibility conditions (but satisfies a suitable generaliza- 
tion). The model employs a piecewise smooth two well potential $(w) with two quadratic 
branches meeting at the spinodal value u = 0. The usual admissibility conditions require 
that the solution vanish (go through the spinodal value) at precisely one point, say £ = 0, 
where it transitions between the two quadratic valleys of so that it must be strictly 
monotone in the neighborhood of the transition point. Instead, the new solutions we find 
are equal to zero on an entire finite interval, and lie in different potential valleys on either 
side of this interval. 

The motivation for considering such solutions comes from a more elaborate model [15] 
where <£>(«) is piecewise quadratic but continuously differentiable. Its graph consists of 
three parabolas, two convex ones separated by a concave one, defined on the spinodal range 
of u values. Any kink solution m(£) that transitions between the two convex potential 
valleys must therefore take values in the spinodal range for £ in some interval, say [—z, z\. 
Here z > is an unknown of the problem and is found to depend on the size of the spinodal 
range. Now $'(«) is trilinear with two increasing branches and a decreasing one between 
them. When the slope of the latter is treated as a parameter and approaches — oo, so 
that the spinodal range tends to degenerate to the point u = 0, a surprising observation 
is made in [15]: The interval — z < £ < z where takes values in the spinodal range 
does not shrink to a point, as one might expect. Rather, z approaches a positive value 
in the limit. The limiting potential is the piecewise smooth biquadratic one considered in 
PQ IU E] where admissible low- velocity kinks have not been found so far. It is thus natural 
to directly try the new type of solutions mentioned in the previous paragraph, that vanish 
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(equal the degenerate spinodal value) over an interval [—z, z], where z > 0, in the context 
of the Atkinson-Cabrera problem. 

While we do not rigorously prove their existence, we are able to construct the new so- 
lutions semi-analytically; we present approximate analytical and also numerical evidence 
of their existence. We find that they bifurcate from the classical Atkinson-Cabrera solu- 
tions, at precisely the threshold velocity Vq below which the latter become inadmissible. 
The new solutions exist essentially over the entire velocity range < V < Vo; thus they 
seem to close the entire gap left by loss of admissibility of the Atkinson-Cabrera solutions. 

The kinetic relation between the applied stress a (constant forcing term) and velocity 
V of moving kinks (dislocations) described by the new solutions is entirely different from 
the (inadmissible) one obtained in [1] over the velocity range < V < Vq] see Fig. [5] 
where the two are compared. Unfortunately, its physical significance is in question, since 
the new solutions appear to be unstable. This is suggested by numerical simulations of 
the Frenkel-Kontorova chain dynamics, where kinks either stop or move with velocities 
above Vq. Recently a similar trend was observed in experiments that measured kinetic 
relations for dislocations in actual two dimensional plasma crystals [10J. Our results are 
consistent with the findings of [15J, where slow kinks are apparently unstable when the 
spinodal range is sufficiently narrow. In contrast, some slow kinks do become stable if 
the spinodal range is sufficiently wide in [15], as is also observed in [5j [TTj . 

The paper is organized as follows. In Section [2] we recall the dynamical driven Frenkel- 
Kontorova model and the equation for traveling wave solutions that describes steady mo- 
tion of a moving dislocation, or kink, under stress. We consider the case of a piecewise 
quadratic potential and recall the main properties of the explicit Atkinson-Cabrera solu- 
tion. The loss of admissibility of this solution in the low-velocity regime motivates us to 
relax the strict monotonicity assumption made in pp. We derive in Section [3] conditions 
for a new type of kink solutions. These conditions include a linear integral equation of the 
first kind, whose kernel is related to the formal Atkinson-Cabrera solution, regardless of 
the admissibility of the latter. Solving the integral equation yields a shape function that 
can be used to obtain the new kink solution. The support [—z, z] of the shape function, 
which depends on the velocity of the moving kink, is the interval where the new kink solu- 
tions take the spinodal value. In Section H] we use linear and quadratic approximations of 
the kernel in the integral equation to obtain approximations of the shape function under 
the assumption that z is small. It turns out that the shape function is a distribution 
that involves two delta functions concentrated at ±z. In contrast, the shape functions 
obtained in [15] are bounded. Our results show that the new solutions bifurcate from the 
Atkinson-Cabrera solutions at the velocity Vq at which the latter become inadmissible. 
In Section [5] we describe the numerical procedure we use to obtain solutions in the case 
when z is not necessarily small. We use this procedure to generate new solutions in the 
low-velocity regime < V < Vq and discuss their properties. We verify the numerical pro- 
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cedure by comparing it with the analytical method described in the Appendix and with 
the analytical results of Section HI In Section [6] we investigate stability of the traveling 
wave solutions using numerical simulations of the discrete chain dynamics. The results 
suggest instability of the new solutions. Section [7] adds viscosity to the model. We find 
that this addition does not seem to stabilize slow new type traveling waves in general. The 
Appendix describes an analytical method for solving a version of the integral equation 
that arises here and in [15], with truncated kernel retaining a finite number of exponential 
terms, by converting it to a differential equation. 

2 Atkinson-Cabrera traveling wave solutions 

The dynamics of the driven Frenkel-Kontorova chain are described by the following equa- 
tion, expressed in dimensionless quantities. 

u n = u n+1 - 2u n + u n -i + n(<j - $'(«„)). (1) 

Here u n (t) is the displacement of the nth mass at time t, u is a ratio of stiffness of 
the nonlinear interaction with the substrate to nearest neighbor interactions and $ is 
the multiple well substrate potential. To model a steadily moving dislocation, we seek 
solutions of (CD) in the form of a traveling wave with (constant) velocity V > 0: 

«n = «(0> i = n-Vt. (2) 
Substituting this ansatz in (DO), we obtain the advance-delay differential equation 

V 2 u" = u(£ + 1) - 2u(0 + u(£ - 1) + u(a - $'(«(£)))• (3) 

We are interested in solutions of that are of kink type. These satisfy the following 
conditions at infinity: 

(«(£)) -»■ u± as £ ±oo, (4) 

where u± are stable constant (uniform) equilibrium solutions of ([3]) located in two different 
wells: 

$'(«±) = cr, U- > u+, &'(u±) > 0. 

The angular brackets in (J3J) denote the average value of the displacement because we 
expect this Hamiltonian discrete system to develop oscillations. The average is taken over 
a sufficiently large interval. 

Instead of the usual periodic potential, we choose a potential with only two wells; 
this is appropriate for describing twinning dislocations. As first shown in pQ, an explicit 



4 



solution of (|3j) , (J4j) can be obtained if one assumes that the substrate potential is piecewise 
quadratic: 

Note that the derivative of this potential is discontinuous at u = 0: 

$'(V) =u — 26 (u) + 1, (6) 

where 9{u) is the unit step function: 9{u) = 1 for u > 0, 0(it) = for u < 0. Observe also 
that in this case one has u± = a ± 1 in OU). 

Suppose that the displacement switches from the second to the first well at £ = 0, so 
that 

u(f) > for £ < 0, < for £ > (7) 

and 

u(0) = 0. (8) 
Then one may replace by #(— £) in (j6]). Then ([3]) becomes a linear equation: 

y V = u{i + 1) - (2 + /i)n(0 + u(£ - 1) + ^(<x - 1 + 20(-£)), (9) 

which can be solved using Fourier transforms. 

It should be emphasized that solutions of fl9]) satisfy the original nonlinear equation 
(j3j if and only if the admissibility condition (jTJ) holds. Otherwise, if a solution of 
violates (JTj) it will be labeled as inadmissible. 

The solution of ([HD constructed by Atkinson and Cabrera [T] (see also [3J, E] for more 
details) is as follows: 



« = u(0 = { 



k£M+(V) Klj k{K, V ) 



(10) 



a + l + 2/i £ ry 77 T/V £<0. 

k£M-(V) Klj k{K, v ) 



Here 

M ± (K) = {& : L(fe, V) = 0, ImJfc ^ 0} U iV^V) (11) 
are the sets of roots of the dispersion relation L(k, V) = 0, where 

L(k,V) =/i + 4sin 2 ^ - V 2 k\ (12) 

with Lk = dL/dk, while 

Ar±(y) = {/c : L(fc, 7) = 0, Im(A;) = 0, kL k (k, V) ^ 0} (13) 
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denote the sets of real roots. The real roots correspond to lattice waves (phonons) emitted 
by the moving dislocation. The construction of N ± (V) implies that corresponding modes 
propagate ahead or behind the dislocation in accordance with the radiation condition 
[U, [9J . This condition, also known as the causality principle [H] , selects solutions such 
that in a frame moving with the dislocation, lattice waves can only be emitted by the 
moving front and must carry energy away from it (thus causing radiative damping). Thus 
phonon modes whose group velocity V g = V + Lk(k,V) / (2Vk) is less than the velocity 
V of the front must be placed behind it (the set N~(V) contributes to £ < 0), while the 
modes with group velocity above V propagate ahead (N + (V)). 

The nonlinearity of the problem is contained in conditions (JTj) and (151) . Applying ([8]), 
one obtains the kinetic relation between the applied stress and the dislocation velocity: 



where the sum is over the set of all real roots, N(V) = N + (V) U N~(V). Thus the applied 
stress is determined entirely by the real roots. As shown in [9J, one can derive (1141) by 
accounting for the energy fluxes carried by the phonon waves. 

For a given V > the solution is thus given by (TTUI) . (TT4"|) . provided that the admissi- 
bility inequalities (j7]) are satisfied. 

Computing the real roots of ffl2|) for each V > 0, we formally obtain the kinetic 
relation ( !T4l) . shown in Fig. UK for the case of fi = 1. The relation consists of disjoint 
segments separated by resonance velocities, i.e. values of V such that L(k, V) = and 
Lk(k, V) = for some real k (see Fig. [lb). A typical solution above the first resonance 
(V = 0.5) is shown in grey in Fig. [2k.. One can see that a moving dislocation emits 
phonon oscillations behind it, with wave number equal to the single positive real root 
of (fl~2l) . As velocity decreases below the first resonance (see the displacement profile at 
V = 0.2 in Fig. [2b), more oscillation modes appear, and w(£) formally obtained from 
(FIU|) . (1141) features phonon emission on both sides. However, closer inspection reveals that 
this solution of ([9]) is in fact inadmissible and should be removed because it violates the 
assumption that u(£) < for £ > (one of the admissibility conditions ([7j)). In fact, 
numerical calculations of the solution U of (Q reveal that all segments of the kinetic 
relation below the first resonance contain only inadmissible traveling waves that change 
signs more than once, and thus need to be removed, while the remaining large-velocity 
segment contains admissible solutions above a certain threshold velocity j9j. In the case 
of fi = 1, solutions of (Q are admissible for V > Vq ~ 0.357. This implies non-existence 
of traveling wave solutions under the assumptions (JTJ), (jU) when the velocity is below 
the threshold value Vq. While this nonexistence issue was recognized by various authors 
[H El El E] , it remains unclear whether there exist any moving kink type solutions below 
the threshold velocity. This paper addresses this question. 




(14) 
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Figure 1: (a) Kinetic relation resulting from the formally obtained Atkinson- Cabrera 
solution. Only the first twelve segments are shown. The grey curves correspond to 
inadmissible traveling waves, (b) Solutions of L(k, V) = for positive real k. The dashed 
lines indicate the first five resonance velocities at which Lk(k, V) = 0. Here fi = 1. 



(a) (b) 




Figure 2: Displacement profiles formally computed from ([TO]) . f|T4|) at (a) V = 0.5 and (b) 
V = 0.2. Solution in (a) satisfies the constraints ((7j) but the one in (b) does not. Here 

/! = 1. 
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3 New traveling wave solutions 



To obtain valid traveling wave solutions of fl3]) in the low-velocity regime, we now replace 
the admissibility conditions ©, (|SJ) by a more general assumption 



u(0 > 0, £ < -z 
u(C) = 0, If I < z 
m(0 < 0, £>z, 



(15) 



where z > is to be determined. In other words, we assume that instead of switching 
wells instanteneously at t — n/V (f = 0), the nth particle may stay at the spinodal value 
u — between the wells over the time period [(n — z)/V, (n + z)/V] (|f| < z). For z = 
we recover Atkinson- Cabrera solutions. 

Adopting a method used in [J5], we observe that $'(«(£)) can now be written as 



$'(u(0) = «(0 + 1 - 2 / h(s)6(s - Z)ds, 

J —z 



(16) 



where we have introduced an unknown shape function h(s), which vanishes outside the 
interval [—z, z] and is normalized so that 



J h(s)ds = 1. 



;i7) 



Thus we obtain 

V 2 u" + (u + 2)u{0 - u{£ + 1) - u(C - 1) = n 



(7-1 + 2 



h{s)0{s - £)ds 



For consistency, we must require that in addition to ( Jig]) and OH), the solution satisfies 
the generalized admissibility conditions f[T5"j) . 

Taking the Fourier transform of (ITS]) and using the convolution theorem, one can show 
(see [7J [16] for details) that 



(19) 



where the kernel is the negative derivative of the solution ([10]) of (j9J) with z = 0: 



?(0 = - 



keM+(V)^k{K, V ) 
ik£, 

keM-(V)^k\K, V ) 



(20) 



8 




At the same time, our assumption that u(£) = at |£| < z implies that = in the 
interval (—z,z). Together with ( 1T91) . this yields the integral equation 

f h(s)q(Z- S )ds = 0, \£\<z. (21) 

Thus the shape function h(£) is an eigenfunction of the integral operator in the left hand 
side of fl2T|) (with kernel given by fl20l) ) associated with the zero eigenvalue. As described 
in more detail in Section HI the integral operator has a zero eigenvalue provided that z 
takes on special values. Note that ( |2T|) is a Fredholm integral equation of the first kind. 
We remark that if instead of ([6]) one considers a trilinear (continuous) $>'(u), 



the same procedure yields a Fredholm integral equation of the second kind, with the right 
hand side of (12T|) replaced by jh(^), where 7 is the width of the spinodal region connecting 
the two wells [TJ)] (see also related problems with different kernels in [7[ [TH] ) . In the limit 
7 — > one recovers (1211) . 

The problem thus reduces to solving the integral equation fl2T]) for z and h(£). Once 
/i(£) and z are known, the convolution theorem yields the traveling wave solution: 

u(0 = 0- £(V) + f h(s)U{£ - s)ds, (22) 

J — z 

where the applied stress 

a = Y>(V)-- f h(s)(U(z-s) + U(-z-s))ds (23) 
2 J—z 

is found by applying u(z) = u(—z) = 0. If z = 0, the shape function reduces to the Dirac 
delta function h(s) = 5(s), and reduce to (ITU]) . pi ]l . respectively. 

We will present numerical evidence that this procedure yields valid traveling waves 
(solutions of (E])) for values of V where the Atkinson- Cabrera solution U of ^ is inad- 
missible. This is due to the fact that for such values of V, m(£) given by (f2"2"|) conforms to 
the generalized admissibility conditions f[T5|) . 



4 Kernel approximations for small z and bifurcation 

Using an approximation of the integral operator in ( 12 ip . we show that the new type of 
traveling waves with z > bifurcate from the Atkinson- Cabrera solutions precisely at 
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the threshold velocity Vo, below which the latter become inadmissible. The conditions 
near bifurcation thus involve values of z near 0. This suggests that in order to study the 
problem analytically, we may replace the kernel q in ( 121 j) on [—z, z] by its piecewise linear 
approximation near zero, as in [7J. We note that g(£) is continuous, while q'(£) has a 
jump discontinuity at £ = 0, as can be shown from (jUJ), since g(£) = —[/'(£). Let 

[qo + q-Z, £ < 0, 

where 

go = g(0), g ± = g'(0±), g + -g_ = 2 / i/^ 2 ; (25) 

the last relation is implied by (Q. Here q is a piecewise linear approximation of q near 
zero. Consider the approximate version of (12~T1) given by J h(s)q(£ — s)ds = 0, |£| < z. 
It turns out that the only L 1 solution of this is the trivial one. To see this, note that 
J_ h(s)q(£ — s)ds can be differentiated twice with respect to £, with second derivative 
equal to (g+ — g_)/i(£) (|^| < z), which therefore has to vanish. In contrast with [15], we 
are led to seek generalized solutions involving delta functions p3]. One way to approach 
this is to study solutions of the corresponding equation of the second kind, which actually 
arises in [T5l. 



( h(s)q(Z-s)ds- 1 h(O=0, |£l<*, (26) 

J — z 

with 7 a positive constant, and then take the limit as 7 — > 0+. Assuming h is smooth 
enough, differentiate ( 1261) twice to find that h must satisfy the ODE 

lh"{® - (g+ - q-)h(0 = 0, |e| < z. 

Inserting the general solution of this, 



h(Z) = Cl eV a + c 2 e-V\ a=V7/(? + -9-), (27) 

into (1261) and evaluating the integral, we obtain an expression linear in £, which must van- 
ish for |^| < z. This yields a homogeneous linear system for the constants C\ and C2, which 
has nontrivial solutions, provided that the corresponding determinant vanishes. This can 
be put into the form of a quadratic equation in the quantity e 2z ^ a , whose coefficients 
depend on z and 7. Solving this quadratic gives the following equations: 

e 2V(9+-'-)/T = B ^ Z > 7 } , (28) 

A{z,i) 

where 

S±(«,7) = -a (qJ + q+ 2 ) 

± a/ (g_ - q+) 2 (g 2 + o 2 (<?- + q+) 2 ) + MoQ-{q- ~ Q+)q+z + 4q_ 2 q + 2 z 2 
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and 

A(z, 7) = qo(q- - q+) + 2q-q+(z - a), 

with a as in ( 127)) . Note that A, B± are continuous algebraic functions, which thus remain 
bounded for finite values of their arguments. This provides a relation between z and 7. 
Suppose that z tends to a positive value as 7 — > 0+. Then the left hand side of ( l28f) . 
hence also the right hand side, grows unbounded. This necessitates that the denominator 
A — > 0, or passing to the limit, A(z, 0) = 0, which reduces to 

qo(q+ ~ q-) 



2q+q- 



(29) 



From (120]) and the last of ([25]) one shows that (q + — g_)/(2g + g_) < 0. Positive solutions z 
of (129)) exist provided —go = U'(0) > 0. On the other hand, the admissibility conditions 
(j7J) for [/ imply that U'(0) < 0. This provides strong analytical evidence that bifurcation 
to the new type of traveling wave with z > occurs precisely at the threshold velocity 
Vo, at which the Atkinson- Cabrera solution becomes inadmissible. In the next section we 
show this numerically as well using the full kernel. 

For small 7 > one can fully determine the constants c\, c 2 and z in terms of 7 by 
enforcing the normalization condition (TTTj) . The result can be put into the form 

m = 9- cosh cosh (*±^ 

a(q- ~ q+) sinh (2£) 

It is straightforward to show that the limit as 7 — > 0+ is 

Kt) = - f ^ M + z)- q+ - 8(£ - z), (30) 

{q- - q+) (9- - q+) 

in the sense of distributions (with delta functions at ±z). One can show this directly by 
substituting the ansatz h(£,) = a+5(£ + z) +a_<5(^ — z) (with a± unknown constants) into 
the first-kind approximate equation f* z h(s)q(^ — s)ds = 0, |£| < z. This determines a± 
and yields (I29|l and (jSP]) . 

One shows that even derivatives of g are continuous at 0. If we add a quadratic 
term g 2 £ 2 to the piecewise linear kernel approximation $M§, the above procedure can be 
repeated with similar results. The analogue of (13"U|) now takes the form 

uic\ uc^ \^ xtt t 2 92 q+-q- -4g 2 2: g+ + g_ 

KO = oc + d(£+z)+a-5(£-z)+(, Q = , a± = —— — =f ; 



Q+-Q-' ' 2 (q+ - q~) 2 (q+ -q- + 4 ^) 

and z is a root of the quartic equation 

(q+ - q-)qo + (4g 2 go - Zq+q-)z — — z - — -z = o. 

3 3(g+ - q_) 
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Note that for g 2 = the last two equations reduces to ( 130)) and (129)) . It turns out that 
only the smallest positive root z of this quartic equation gives rise to traveling waves 
satisfying the generalized admissibility conditions ( TH>|) ; this root is well approximated by 
that of (129)). 

Another possibility is to consider a kernel in the form of a sum of exponentials, but 
containing only a finite number N of terms from ( |20l) . In that case, ( l2Tj) can be solved 
analytically. The procedure is described in detail in the Appendix and is not confined 
to the case of small z. Unfortunately, this approach involves matrix computations that 
become progressively more difficult computationally as N is increased. Instead, in the 
following section we solve (|2T|) numerically, using the analytical method of the Appendix 
to validate our numerical computations. 

5 Numerical results 

We now consider the general case where z is not necessarily small. 

For given V > 0, the solution and z of ()2"Tj) was found numerically as follows. Since 
there are infinitely many roots of (TT2")) . it was necessary to approximate the kernel (I20p by 
retaining the first iV roots of ffl2l nearest the origin. To get an accurate approximation of 
the kernel g(£) near £ = 0, it was necessary to include a large number of roots (N = 400 
was typically used). The trapezoidal approximation of the integral equation for a finite z 
(with 400 uniformly distributed mesh points) then yielded a homogeneous linear system 
Q(z)h = 0. Here h is the vector of values of h(s) at the mesh points. Solving detQ(^) = 
for z, we obtained the corresponding h normalized to satisfy ( ITT)) in the sense of the 
trapezoidal approximation. In general, there were more than one root of detQ(,z) = 0, but 
at most one of these yielded admissible solutions that satisfied the generalized admissibility 
conditions (I15p within numerical error. Once the admissible z and h were found, the 
trapezoidal approximation of the integrals in ( 12 2 j) and ( 12 3 j) was used to compute the 
solution and the applied stress a. Fig. [3] shows the solution it(£) obtained at V = 0.2 
and n — 1 (black curve) along with the inadmissible Atkinson- Cabrera solution £/(£) 
(z = 0, grey curve). In this case z = 0.21. One can see that the obtained solution satisfies 
the constraints (TT5)) within numerical error (which is of the order of 10~ 7 in this case). 

Repeating this procedure for a range of velocities, we obtain the kinetic relation o~(V) 
and the corresponding relation z(V) between z and V shown in Fig. |H Note that z(V) 
decreases as V grows and becomes zero at V = Vq ~ 0.357. This is the threshold velocity 
such that g(0) = 0, which equals the bifurcation velocity as predicted in Sec. HJ At 
V > Vq, the Atkinson- Cabrera solutions are admissible, and we have z = and a = S(V) 
(thick segments in Fig. H]). Interestingly, for V below Vq we obtained admissible solutions 
with z > even in the immediate vicinity of the resonance velocities. At these velocities, 
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(a) (b) 




Figure 3: (a) Traveling wave solution at V = 0.2, fi = 1 with z = 0.21 (black curve), 
shown together with inadmissible z = solution U (grey curve), (b) Zoom-in of the 
z = 0.21 solution inside the rectangle in (a). 



(a) (b) 
AC solutions (z=0) 7 




Figure 4: (a) Kinetics relation cr(V) and (b) the corresponding z(V) at /i = 1. The thicker 
segments above the threshold velocity Vo ~ 0.357 indicate the parts of the curves that 
correspond to Atkinson- Cabrera solutions (z = 0). 



13 




Ol 0.2 03 OA 0.5 y 



Figure 5: Comparison of the kinetic relation a(V) (black curve) generated by the new 
solutions to the relation S(V) (grey lines) obtained from Atkinson- Cabrera solutions. The 
two curves coincide above the threshold velocity Vo ~ 0.357. The grey curve corresponds 
to inadmissible Atkinson- Cabrera solutions below Vq. Here \l — 1. 

the Atkinson-Cabrera kinetic relation S(V) is inadmissible and has singularities, while 
the kinetic relation <j(V) based on the present admissible solutions has cusps; see Fig. [5j 

The computed shape function h(s) corresponding to velocity V = 0.353 (z = 0.0026) 
is shown in Fig. EK Note that it approximates the delta-function behavior at s = ±z 
that was predicted in Sec. HI As V decreases, the shape functions develop oscillations due 
to the increasingly oscillatory behavior of g(£) at smaller V; see, for example, Fig. [6b for 
h(s) at V = 0.2, which corresponds to the solution shown in Fig. [3j 

To see the effect of the number of roots included in the truncated kernel, we re- 
peated the simulation with N = 40 roots. The resulting kinetic relations are compared 
in Fig. [Tk^b. One can see that the two relations are very close everywhere except at 
quite small velocities V < 0.05, where inclusion of only 40 roots does not approximate 
the Atkinson-Cabrera solution well enough. In contrast, there is a noticeable difference 
between the corresponding z(V) curves (see Fig. [Tfc). The lower values of z when fewer 
roots are included are due to the slow convergence of the kernel g(£) in (12"U|) near £ = 0. In 
particular, one obtains a lower threshold value Vq, where z becomes zero, as the number 
of roots in the truncated kernel is decreased. 

Nevertheless, we can use computations with a small number of roots to verify our 
numerical procedure by comparing shape functions computed numerically with the ones 
obtained using the semi-analytical method described in the Appendix. The results for 
N = 8 and iV = 44 roots are shown in Fig. [SJ where the numerical results (thin curves) 
are compared to the regular part h (x) of the semi-analytical shape function h(x) = 
a + 5(x + z) + a-5(x — z) + ho(x), shown by the thick curve, at V = 0.245. One can see 
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Figure 6: Shape functions h(s) in the interval [—z, z] at (a) V = 0.353, z = 0.0026 and 
(b) V = 0.2, z = 0.21. Here y.= l. 

that the shape functions obtained using the two methods agree very well (the difference 
is below 0.0004% in (a) and 0.007% in (b)). Observe, however, that including more roots 
resulted in a more oscillatory function and, as remarked above, higher value of z. 

Finally, we consider the bifurcation diagram for z(V) near the threshold velocity Vq. 
To obtain accurate results near this threshold, one needs to include an even larger number 
of roots in the truncated kernel that approximates ( 1201) . Instead, we compute g(£) using 
numerical evaluation of the integral representation of the kernel without truncating roots, 
and compare the resulting zj^(V) obtained using the trapezoidal rule (solid curve) to the 
curve zl(V) obtained from (125]) using the linear approximation of the kernel (dashed line) 
in Fig. [9j As one can see, the two curves become closer as we approach the threshold 
velocity from below. 

6 Stability of traveling wave solutions 

To investigate the stability of the new type of traveling waves, we conducted a series of 
numerical simulations, since in general it is quite difficult to check stability of traveling 
waves analytically [2]. For a given applied stress a, we used the Verlet algorithm (a 
symplectic scheme) to solve the system (TjQ) of ordinary differential equation for a truncated 
lattice with N masses, ranging from N = 600 to 2000, depending on the time of the 
simulation. A longer chain was used if the simulation ran for a long time, in order to 
avoid reflection of elastic waves from the domain boundaries. The boundary conditions 
were uq = cr + 1 and un = o — 1. Two types of initial conditions were used. The first one 
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(c) 



Figure 7: (a) Comparison of the numerical results at \x = 1 obtained using the truncated 
kernel with N = 400 (solid curves) and N = 40 roots (dashed curves): (a) kinetic relations 
<j(V); (b) the difference between the two relations; (c) the corresponding z(V) curves. 
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Figure 8: Comparison of the numerically computed shape functions h(s) in the interval 
[—z, z] at V = 0.245 (thin curves) and the regular part of the shape functions obtained 
using the semi- analytical method described in the Appendix (thick curves) including (a) 
8 roots (z = 0.1006) and (b) 44 roots (z = 0.1363). Here /i = 1. 



Z Zl~Zn 




0.345 0.350 0.355 V 



(a) (b) 

Figure 9: Comparison of the function z^{V) computed numerically using the full kernel 
(solid line) and zl(V) obtained from (125]) using the linear approximation of the kernel 
(dashed line) near the bifurcation point Vq ~ 0.357. Here // = 1. 
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was Riemann-type piecewise constant initial displacement and zero initial velocity: 

{(a + 1,0) < n < n 
(0,0), n = n , (31) 

(a -1,0), n <n<N 

where n = N/2 for even N. Numerical simulations with these initial data sought to 
identify stable states at a given loading that have a relatively wide basin of attraction. To 
capture other possibly stable states that coexist with solutions found using (13T1) but have 
a more narrow basin of attraction, and to identify solutions that are likely to be unstable, 
we used a second type of initial conditions, that were built from the obtained traveling 
wave profiles u n {t) — u(n — Vt): 

{(a + 1,0) 0<n<p 
(u(n-n ),-Vu'(n-n )) p < n < N - p. (32) 
(a -1,0) N-p<n<N 

The truncated traveling wave solutions were surrounded by intervals of constant dis- 
placement of appropriately chosen size p < n in order to ensure compatibility with the 
boundary conditions and avoid wave reflection from the boundaries. In both types of sim- 
ulations, after a sufficiently long time, the solution approached an attractor corresponding 
to either a stationary dislocation (zero velocity) or a steadily moving front. 

The results are shown in Fig. [TOj One can see that when the applied stress is below a 
certain threshold <7d (here o"d = 0.128), the long-time solution features a stationary front. 
For example, in the simulations with piecewise-constant initial conditions (l3Tj) the front 
propagates for some time (which increases as we approach the <jd from below) and then 
becomes stationary. This is ullustrated in Fig. [HI which shows the position u(t) of the 
front (defined as the value of n such that u n {t) > and u n+ i(t) < 0) at a = 0.125 and 
cr = 0.1275. In general, stable stationary solutions exist when a is inside the trapping 
region \a\ < a P , where a P = [i/(A + /i) (^ 0.447 for ji — 1) is the Peierls stress [P7] . 
The fact that o^, < op has been observed in earlier works, e.g. [5j [3]. The trapping region 
is marked by a thick segment along V = in Fig. [TUl 

When stress is above the threshold value (a > a D ), the solution approaches a steady 
dislocation motion after some time. For example, at a — 0.14, the long-time solution 
features the dislocation moving steadily with velocity V = 0.54; see Fig. [12j Comparison 
of the numerical solution zoomed in around the front (circles) and the corresponding 
traveling wave solution (solid curve) in Fig. [T2b shows excellent agreement, indicating 
stability of the traveling wave. In general, our simulations indicate that at stresses above 
o"d all traveling waves are stable. Note that the threshold value <td (the dynamic Peirels 
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Figure 10: (a) Results of the numerical simulations at fi — 1 with initial data fl3Tl) (black 
circles) and (1321) (grey circles), shown together with the kinetic curve, (b) Zoom- in inside 
the rectangle in part (a). Thick black segment along V = axis indicates the trapping 
region. Thinner portion of the kinetic curve corresponds to the solutions with z > 0. 



(a) (b) 




Figure 11: Position v(t) of the dislocation at (a) a = 0.125 and (b) a = 0.1275 in the 
numerical simulations. Here \i — \. 
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(a) (b) 




Figure 12: (a) Displacement profiles (solid lines) at t — and t = 400 for the numerical 
simulation at \i — 1 and a = 0.14. (b) The numerical solution (circles) at t = 600 zoomed 
in around the dislocation front and compared to the traveling wave solution (solid curve) 
with the same velocity, V = 0.54. 

stress) corresponds to a local minimum of the kinetic curve, as hypothesized in [T], and 
is below the Peierls stress dp, implying that stable stationary states and stable steady 
motion coexist at stresses between the two values. Similar stability results were reported 
in [3] for simulations that include a viscosity term. A proof of stability of traveling waves 
with sufficiently high velocities, which for technical reasons does not extend to the whole 
o~ > o~d region, can be found in [2J. 

Note that inside the stability interval a > ctd suggested by the numerical simulations, 
all traveling waves have z — 0. Under initial data fl32|) based on the new-type traveling 
wave solutions with z > 0, numerical simulations always converged to attractors with 
stationary fronts. Our results thus suggest that such solutions are likely to be unstable. 
We remark, however, that some low-velocity solutions apparently become stable when a 
sufficiently wide spinodal region is included [TT| [T5]. 

7 Effect of viscosity 

We now consider the effect of adding viscosity to the model on kink solutions and their 
stability. The rescaled governing equations become 

u n + au n = u n+ i - 2u n + w n __i + fi(a - $'(«„)), (33) 

where a > is the dimensionless viscosity coefficient. Traveling wave solutions are given 
by (flOjh (PI) for z = and by (122]), p5]) for z > 0, with flTJJ replaced by 

L(k, V) = /X + 4 sin 2 - - V 2 k 2 - ikaV. (34) 
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Figure 13: (a) Kinetic relation a(V) at ji — 
ficient a. (b) Zoom-in on the curves with a 
kinetic curves correspond to solutions with z 
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and different values of the viscosity coef- 
and a = 0.01. Thicker portions of the 
0. 



In this case there are no real roots; at small a the roots shift away from the real axis into 
upper and lower halves of the complex plane according to the radiation condition [9]. 

The effect of viscosity on the kinetic relation is shown in Fig. [13j As expected, viscosity 
smoothens the cusps at the resonance velocities. The amplitude of oscillations that are 
very pronounced in the kinetic curves at small a decreases as a increases; the velocity V 
at which z = solutions become admissible decreases as well. See also the corresponding 
z(V) graph in Fig. [T4"r . At sufficiently large a we have Vq = 0, so that all z — solutions 
become admissible. 

The value of Vq also becomes smaller as we decrease /z, as shown in Fig. [14b . Another 
effect of the parameter /i is the different values of the resonance velocities. 

The stability picture is not significantly affected by a and \i. See, for example, Fig. [T5l 
where /i = 0.5 and a = 0.1. As before, only sufficiently fast z = solutions, above the 
last local minimum of the kinetic curve, appear to be stable. 
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Figure 14: (a) z(V) at fi — 1 and different values of the viscosity coefficient a corre- 
sponding to the kinetic curves in Fig. [T3J (b) z(V) at a = 0.1 and different values of 




Figure 15: Results of the numerical simulations at [i = 0.5 and a = 0.1 with initial data 
( I3T1) (black circles) and ( 132|) (grey circles), shown together with the kinetic curve. Thinner 
portion of the kinetic curve corresponds to the solutions with z > 0. 
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A Appendix: Exact solution of the integral equation 
with truncated exponential kernel 

Suppose that in ffTTj) we truncate the roots so that we keep an equal number of complex 
roots in M + (V) and M~(V). Due to the presence of real roots, M + (V) and M~(V) will 
contain a different number of roots of L(k, V) = 0. Suppose that after truncation we keep 
a total number N of roots K n , M of which are in M + (V) and N — M are in M~(V), i.e., 

JM+(Y) forn = l,...,M, 
n [M-(V) for n = M + l,...,N 

Let 

A n = 2^/L K {n n ,V) (A.2) 



and define, in view of (120)) . 

?+(£) = E ^e fc "«, £ > 0, 

9(0 = { n= N (A.3) 

?-(0= E (-^n)e fc "«, £<0. 

n=M+l 

Because is not smooth at 0, we observe that 

g(£ — s)h(s)ds — / — s)h(s)ds + / q_(£ — s)h(s)ds, 



j-z Ji 

and each of the last two integrals has a C°° kernel and can be differentiated. The same is 
true for the integrals 

J n (0 = / A n e k ^h(s)ds, n = 1, . . . , M, 

J - r Z z (A.4) 
J n (0 = / (- A n )e k ^-^h(s)ds, n = M + 1, . . . , N. 



In view of them and (1A.2|) . (1A. 3[) . the integral equation f* q(£ — s)h(s)ds = /(£) becomes 



N 



= /(O, ki<*, (a.5) 



n=l 



where / is either identically zero or a given function or 7/1 for the problem considered in 
[T5] and here in Section HI Differentiate (1A.4I) with respect to £ using the Leibnitz Rule 
to find 

r n (0=A n h(0 + k n I n (0. (A.6) 
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Then differentiate the integral equation (I A. 5 1) and use ( 1A.6I) to eliminate I' n . Repeat 
N — 1 times so that together with (I A. 5 1) we have A" equations where we have used ( 1A.6I) 
to eliminate I' n after each differentiation |12| . Letting 

/ (m) (0 = ^ m /(0M m 

for m = 0, 1, . . . , the resulting system of A" equations for I n is 



N 



m—1 



EcX(o = / K1) (o-E 

P =i 



n=l 



N 



n=l 



m 



N. 



(A.7) 



For m = 1 the right sum above is understood to vanish. The right-hand side is a linear 
combination of derivatives of h and /. The coefficient matrix of /„ in the left-hand side 
is the N x N Vandermonde matrix evaluated at the roots: 





1 


l 


1 




h 




k N 


[Vmn] = K- 1 } = 


k\ 


k\ 


k 2 




k?- 1 


h.N-1 . 

hv 2 


■ k N ~ l 



(A.8) 



Whenever all roots are distinct this matrix is invertible. Solving ( 1A.7j) for /„ yields 
expressions of the form 



I n = L n [h]-L n [f], n = l,...,N, 



(A.9) 



where L n and L n (n = 1, . . . , N) are linear differential operators of order N — 2 and N — 1, 
respectively, with constant coefficients depending on k m and A m . Differentiate the last of 
flA.7j) (m = N) once more and use flA.61) to eliminate I' n to find 



N N 

p=i L 



n=l 



N 



n=l 



(A.10) 



Substitute (1A.9|) into flA.lOp . collect terms with h and its derivatives on the left, / and 
its derivatives on the right to arrive at an equation of the form 



L[h]=G[f], 



(A.11) 



where L is a linear differential operator of order at most A^ — 1 and G is a linear differential 
operator of order N, both with constant coefficients depending on k m and A m (shown 
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explicitly in ( 1A.23I) and (I A. 241) below). Thus any solution h of ( 1A.5I) satisfies the ODE 
(lA.lip of order N — 1. Next note from flA.4jl that 



l n (-z)=0, n = l,...,M, 
Uz)=0, n = M + 1,...,N. v ; 



Together with flA.91) this yields A" boundary conditions for h: 

L n [h](-z) =L n \f](-z), n=l,...,M, 
L n [/i](*) =L n [/M n = Af+l,...,iV. 

These are of order N — 2 for h (shown explicitly in (IA.25j) below). Hence any solution h of 
flA.5j) necessarily satisfies ODE flA.llj) with boundary conditions (IA.13I) . Conversely, let 
h satisfy flATT]) . flA~T3|) . Define 4 by flA~9|) . Then they satisfy flA~7l) . flA~T0l) and (jXl2j> . 



Differentiate flA.7[) and subtract from it (IA.TD with m replaced by m + 1 (subtract (IA.10I) 
if m = N) to find after some calculation 



N 



n=l 

This is a homogeneous linear system for the terms in brackets, whose coefficient matrix 
is the invertible Vandermonde matrix. Hence the term in brackets above must vanish for 
each m, therefore (1A. 6|) holds. The unique solution of the I VP flA.6j) . flA.12j) for I n (£) 



with h given is ( 1A.4I) . But since /„, satisfies ( 1A.9I) , setting m = 1 in the latter yields ( 1A.5I) . 
VKe Ziave shown that h G C Af_1 (— z, z) is a solution of ( 1A.5I) z/ and on/y «/ it satisfies the 
2-pomtBVP (IA~TT|) . flA~T3l) . 



Observe that the ODE is of order at most N — 1 for /i, but there are A^ boundary 
conditions. Consider the eigenvalue problem of replacing /(£) with 7/1 (^) in (1A.5j) . (1A.11I) . 
(1A. 13h . Recalling that G is an N th order operator, we obtain an A^ th order ODE for h, 
L[h] — ^G[h] = with A^ boundary conditions of order N — 1 and the difficulty is removed. 
Note that in this ODE 7 multiplies the highest (A^ th ) derivative of h, so the limit as 7 — > 
is to be taken with care. The A^ constants of the general solution of L[h] — ^G[h] = now 
satisfy N homogeneous linear equations (from flA.131) ) which only have the trivial solution 
unless the determinant of the coefficient matrix vanishes, which is an algebraic equation of 
the form ^(7, z) = (the constants are then nonunique unless the normalization condition 
ffT7|) is enforced.) 

Presumably '^(7, z) = can be solved for 7 in terms of z. We showed in Section |4] for 
a piecewise linear kernel, that as 7 — > and z approaches a positive value zo, such that 
i/j(0, Zq) = 0, h, which is a combination of real exponentials, tends to two delta functions 
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at ±z. Assuming that this also occurs here, suppose that the solution to the homogeneous 
integral equation 

N 

^4(0 = 0, \£\<z (A.14) 

n=l 

is of the form 

h(s) = a + 5(s + z) + a-5(s - z) + h (s), |f| < z, (A.15) 

where 5 is a delta function and h is suffieciently smooth. For the piecewise quadratic 
kernel considered in Section H] , h Q is a constant. Here a± are unknown constants while 
the normalization condition (I17j) becomes 

/h (s)ds = 1 — a + — ct_. (A. 16) 

-z 

In view of (1A.15j) ho in place of h in ( 1A.5j) solves the integral equation 

N M N 

E J «(0 = "«+ E 4»e**« + *> -a- E (-4.)e*"«-\ |£| < z. (A.17) 

n=l n=l n=M+l 

Suppose Sj, = 0, . . . , N, is the j th elementary symmetric function of k — (hi, . . . , k n ) e 

iV 

R N (So = 1, 5*1 = X] kn> ••• > Sjv = rin=i^n)- Namely, Sj are the fundamental scalar 

n=l 

invariants of the diagonal N x N matrix diag(k) with diagonal entries fcj. Then its 
characteristic polynomial is 

N N 

P(x) = det (xl - diag(k)) = Y[(x-k n ) = ^2(-l) N - n S N -nX n , xeTR, (A.18) 

n=l n=0 

where I is the N x N identity matrix. Let K = N — 1 and 

= (^1; • • • i ki—l, fci-fl, • • • , fc n ) £ R- 

(obtained by removing ki from k). The characteristic polynomial of the K x K matrix 
diag(k /; ) is 

^ X l > n=0 

Here S(j^ is the j th elementary symmetric function of kn with £(0,*) = 1 and I is the 
K x K identity matrix. Then the inverse of the Vandermonde matrix (4) has components 

- ^wtr - (A - i9) 
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In particular, letting the i th Lagrange polynomial be k(x) = P{(x) / Pi(ki) , Wij equals the 
coefficient of x^ 1 in li{x): 



N 



l i (x) = ^2w ij x i - 1 



(A.20) 



It can be shown that 



N UN o 

EK i 0(7V-m,i) _ c 



. m , m = l,...,N. 



A' 



From this and (1A.19I) it follows that the solution J n of ^2 k™ 1 J n = d m satisfies 



n=l 



N 



N 



+1— m""m- 



(A.21) 



n=l 



m=l 



Next we write the differential operators in (lA.llj) . ( 1A.13I) more explicitly. Letting d r , 
equal the right-hand side of (1A.7I) . and letting 



N 



v = 0, v 9 = £ k t lA ni q=l,...N, 



(A.22) 



n=l 



substitute into (IA.21I) with m = p + 1 and subtract the result from flA.101) to find 



N 



£(-1)"-^ 

p=0 



p-l 



n=0 



(A.23) 



or 



N 



p-i 



N 



(A.24) 



p=0 n=0 p=0 

The above is the explicit form of the differential equation (lA.llj) . Observe that it is of 
order N — 1 in h but order N in /. 

Next we express the boundary conditions explicitly. Recall that the solution Jj of 



J2n=i k } X Jj = d% is Ji = J2f=iWijdj with as in (17). The boundary conditions 
(TAT21) become, with K = N - 1, 



A 



£(-l) X P S(K-p,i) 
p=0 
K 

?X~ l ) K ~ PS {K-p,i) 



p=0 



p-l 

Y,Vp-*h in H-z)-fto( 

,n=0 
p-l 



x;«^»a (b) («)-/ 



(p) ( 



n=0 



= 0, z = l,...,M, 
0, i = M + l,...,A. 



(A.25) 
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We now consider the symmetries of the root sets M ± (V), N ± (V) in (fTTT) and ( 1131) . 
Let Mfiy) = M ± (V) - N±(y). Then (recall fc„ = m n ) k e M+(V) if and only if 
—k G M~(V). Also k G M± (V) if and only if -k G M±(V), respectively. In addition 
k G N ± (V) if and only if — k G N ± (V), respectively. Suppose we keep these symmetries in 
the truncated set of roots. Hence for some integers R±, C we will have: M = 2R + +2C + 1 
is the number of roots in M + (V), N — M = 2i?_ + 2C+l is the number of roots in M~(V). 
Here R± is the number of real positive roots in N ± (V), respectively; C is the number 
of complex roots in the first quadrant and there is one purely imaginary root in each of 
M ± {V); see [9] for details. Recall that q±(£) and their even derivatives are equal at £ = 0. 
This implies that all odd coefficients v m in (1A.22j) vanish: 



N 



v 2q+1 =Y,k 2 n q A n = 0, g = 0,l,... (A.26) 



n=l 



In particular, Yln=i = and the coefficient of the highest order derivative h£ N ^ in 
flA"^3l . flA~24l) and h^ N ~^ in flA~25|) vanishes. Hence ODE f[Q4D is of order N — 2 va. h 
and BC (IA.25[) of order — 3 (but still of order A^ and AT — 1, respectively, in /). This 
is important. The general solution of the ODE will have N — 2 undetermined constants, 
which together with a± in (1A.15I) and z comprise A^ + 1 constants. These are to satisfy 
the N + 1 conditions (EP5i) . (jATTBIl . 

We turn to the problem for the regular part ho of h; see (1A.15I) through (1A.17I) . Letting 
D n f = f n \ the right hand side of (lATTTjl or flAT24l) above is G[f] = P(D)[f] where P 
is the polynomial in (1A. 181) .Observe that for / equal to the right-hand side of ( 1A.17I) . 
G[f] = since P(k n ) = 0. Hence flA.24j) becomes homogeneous: 



N p-2 

J2(-^ N - p S N -pJ2 v P~nh { o\0 = 0, -z<i<z. (A.27) 

p=0 n=0 

The boundary conditions do not; however in the i th equation of ()A.25j) the term involving 
/ is Pi(D)[f] (see (jX|) : here D is the derivative operator). Note that Pi(kj) = if i ^ j, 
hence we obtain (with K = N — 1): 

K p-2 

Y,h~L) K ~ P S(K- P ,i) J2 v P-nh { n \-z) = -a+AiPifc), i = 1, . . . , M, 

p=0 n=0 

K p-2 

Y,(-V K ~ P S { K-p, t ) v P - n h { n Xz) = -a-i-A^Pih), i = M + 1, . . . , N. 

p=0 n=0 



(A.28) 



In (1A.27I) . ( 1A.28I) the summation over n has upper limit p — 2 in view of (1A.26I) . It is 
understood that the sum is zero for p < 2. 
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The N conditions (1A.28I) comprise a homogeneous system for the N — 2 coefficients of 
the general solution ho of ( 1A.27I) together with a ± . For nontrivial solutions a determinant 
involving z has to vanish; this determines z. The resulting loss of uniqueness in the 
solution of the system is hopefully removed by the additional equation flA.16j) . 
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